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^vq Abstract 

Quantization of compressed sensing measurements is typically justified by the ro- 
bust recovery results of Candes, Romberg and Tao, and of Donoho. These results 
guarantee that if a uniform quantizer of step size 5 is used to quantize m measure- 
ments y = <3?x of a fe-sparse signal x € M. N , where satisfies the restricted isometry 
property, then the approximate recovery via ^-minimization is within 0{5) of x. 
The simplest and commonly assumed approach is to quantize each measurement inde- 
c/3 pendently. In this paper, we show that if instead an rth order SA quantization scheme 

i with the same output alphabet is used to quantize y, then there is an alternative recov- 

ery method via Sobolev dual frames which guarantees a reduction of the approximation 
error by a factor of (ra/fc)^ -1 / 2 )" for any < a < 1, if m > r A;(log iV) 1 ^ 1- "). The 

(N) result holds with high probability on the initial draw of the measurement matrix $ 

from the Gaussian distribution, and uniformly for all /c-sparse signals x that satisfy a 

^ mild size condition on their supports. 

1 Introduction 

Compressed sensing is concerned with when and how sparse signals can be recovered exactly 
or approximately from few linear measurements (8| JTo|[l4| . Let $ be an m x iV matrix 
providing the measurements where m N, and Eydenote the space of /c-sparse signals 
in R , k < m. A standard objective, after a suitable change of basis, is that the mapping 
x i — y y — $x be injective on . Minimal conditions on $ that offer such a guarantee 
are well-known (see, e.g. Ill]) and require at least that m > 2k. On the other hand, 
under stricter conditions on $, such as the restricted isometry property (RIP), one can 
recover sparse vectors from their measurements by numerically efficient methods, such as 
^-minimization. Moreover, the recovery will also be robust when the measurements are 
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corrupted [9], cf. 15 ; if y = $x + e where e is any vector such that ||e|| 2 < e, then the 
solution x^of the optimization problem 

min ||z||i subject to \\$z — y\\2 < e (1) 

will satisfy ||sc — x # ||2 < e. 

The price paid for these stronger recovery guarantees is the somewhat smaller range of 
values available for the dimensional parameters m, k, and N. While there are some explicit 
(deterministic) constructions of measurement matrices with stable recovery guarantees, best 
results (widest range of values) have been found via random families of matrices. For exam- 
ple, if the entries of $ are independently sampled from the Gaussian distribution A/"(0, — ), 
then with high probability, $ will satisfy the RIP (with a suitable set of parameters) if 
m ~ fclog(y). Significant effort has been put on understanding the phase transition behav- 
ior of the RIP parameters for other random families, e.g., Bernoulli matrices and random 
Fourier samplers. 



Quantization for compressed sensing measurements 

The robust recovery result mentioned above is essential to the practicality of compressed 
sensing, especially from an analog-to-digital conversion point of view. If a discrete alphabet 
A, such as A = SZ for some step size S > 0, is to be employed to replace each measurement 
Uj with a quantized measurement qj :— j)j G A, then the temptation, in light of this result, 
would be to minimize ||e||2 = \\y — q\\z over q G A m . This immediately reduces to minimizing 
\yj — qj\ for each j, i.e., quantizing each measurement separately to the nearest element of 
A, which is usually called Pulse Code Modulation (PCM). 

Since \\y — q\\ 2 < \8y/m, the robust recovery result guarantees that 

\\ x ~ x *cuh ^ <V™- (2) 

Note that ^ is somewhat surprising as the reconstruction error bound does not improve 
by increasing the number of (quantized) measurements; on the contrary, it deteriorates. 
However, the y/m term is an artifact of our choice of normalization for the measurement 
matrix $. In the compressed sensing literature, it is conventional to normalize a (random) 
measurement matrix $ so that it has unit-norm columns (in expectation). This is the 
necessary scaling to achieve isometry, and for random matrices it ensures that E||$x|| 2 = ||x|| 2 
for any x, which then leads to the RIP through concentration of measure and finally to the 
robust recovery result stated in ([T]). On the other hand, this normalization imposes an in- 
dependent dynamic range for the measurements which scales as 1/ \/m, hence it is not fair to 
use the same value 5 for the quantizer resolution as m increases. In this paper, we investigate 
the dependence of the recovery error on the number of quantized measurements where 5 is 
independent of m. A fair assessment of this dependence can be made only if the dynamic 
range of each measurement is kept constant while increasing the number of measurements. 
This suggests that the natural normalization in our setting should ensure that the entries of 
the measurement matrix $ are independent of m. In the specific case of random matrices, we 
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can achieve this by choosing the entries of $ standard i.i.d. random variables, e.g. according 
to jV(0, 1). With this normalization of $, the robust recovery result of [9], given above, can 
be modified as 

\\y-yh<e => \\x - x*\ 2 < -7=e, (3) 

i/m 

which also replaces (|2]) with 

Ik-^pcMlh £ 5- (4) 

As expected, this error bound does not deteriorate with m anymore. In this paper, we 
will adopt this normalization convention and work with the standard Gaussian distribution 
A/"(0, 1) when quantization is involved, but also use the more typical normalization Af(0, 1/m) 
for certain concentration estimates that will be derived in Section The transition between 
these two conventions is of course trivial. 

The above analysis of quantization error is based on PCM, which involves separate (inde- 
pendent) quantization of each measurement. The vast logarithmic reduction of the ambient 
dimension N would seem to suggest that this strategy is essentially optimal since information 
appears to be squeezed (compressed) into few uncorrelated measurements. Perhaps for this 
reason, the existing literature on quantization of compressed sensing measurements focused 
mainly on alternative reconstruction methods from PCM-quantized measurements and vari- 
ants thereof, e.g., j§@0|2^[22|[29]. The only exception we are aware of is u\, which uses 
EA modulation to quantize x before the random measurements are made. 

On the other hand, it is clear that if (once) the support of the signal is known (recovered), 
then the m measurements that have been taken are highly redundant compared to the max- 
imum k degrees of freedom that the signal has on its support. At this point, the signal may 
be considered oversampled. However, the error bound Q does not offer an improvement of 
reconstruction accuracy, even if additional samples become available. (The RIP parameters 
of $ are likely to improve as m increases, but this does not seem to reflect on the implicit 
constant factor in Q satisfactorily.) This is contrary to the conventional wisdom in the the- 
ory and practice of oversampled quantization in A/D conversion where reconstruction error 
decreases as the sampling rate increases, especially with the use of quantization algorithms 
specially geared for the reconstruction procedure. The main goal of this paper is to show 
how this can be done in the compressed sensing setting as well. 



Quantization for oversampled data 

Methods of quantization have long been studied for oversampled data conversion. Sigma- 
delta (EA) quantization (modulation), for instance, is the dominant method of A/D con- 



version for audio signals and relies heavily on oversampling, see 13,18, 24] . In this setting, 
oversampling is typically exploited to employ very coarse quantization (e.g., 1 bit/sample), 
however, the working principle of EA quantization is applicable to any quantization alpha- 
bet. In fact, it is more natural to consider EA quantization as a "noisa^] shaping" method, 

a The quantization error is often modeled as white noise in signal processing, hence the terminology. 
However our treatment of quantization error in this paper is entirely deterministic. 
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for it seeks a quantized signal (qj) by a recursive procedure to push the quantization error 
signal y — q towards an unoccupied portion of the signal spectrum. In the case of bandlimited 
signals, this would correspond to high frequency bands. 

As the canonical example, the standard first-order EA quantizer computes a bounded 
solution (uj) to the difference equation 

(Au)j := Uj-u j - 1 = yj-qj. (5) 

This can be achieved recursively by choosing, for example, 

qj = argmin + Uj — p\. (6) 

Since the reconstruction of oversampled bandlimited signals can be achieved with a low- 
pass filter ip that can also be arranged to be well-localized in time, the reconstruction error 
ip * (y — q) = Aip * u becomes small due to the smoothness of <p. It turns out that, with 
this procedure, the reconstruction error is reduced by a factor of the oversampling ratio A, 
defined to be the ratio of the actual sampling rate to the bandwidth of <p. 

This principle can be iterated to set up higher-order EA quantization schemes. It is well- 
known that a reconstruction accuracy of order 0(A _r ) can be achieved (in the supremum 
norm) if a bounded solution to the equation A r u = y — q can be found [13] (here, r G N is 
the order of the associated EA scheme). The boundedness of u is important for practical 
implementation, but it is also important for the error bound. The implicit constant in this 
bound depends on r as well as ||w||oo- Fine analyses of carefully designed schemes have shown 
that optimizing the order can even yield exponential accuracy 0(e~ cA ) for fixed sized finite 
alphabets A (see (l8j), which is optimal apart from the value of the constant c. For infinite 
alphabets, there is no theoretical lower bound for the quantization error as A increases. 
(However almost all practical coding schemes use some form of finite alphabet.) 

The above formulation of noise-shaping for oversampled data conversion generalizes nat- 
urally to the problem of quantization of arbitrary frame expansions, e.g., (3). Specifically, 
we will consider finite frames in IR fc . Let E be a full- rank m x k matrix and F be any left 
inverse of E. In frame theory, one refers to the collection of the rows of E as the analysis 
frame and the columns of F as the synthesis (dual) frame. For any x G M fc , let y = Ex be its 
frame coefficient vector, q e A m be its quantization, and let x := Fq be its reconstruction 
using the dual frame. Typically A m C\y + Ker(F) = 0, so we have x ^ x. The reconstruction 
error is given by 

x-x = F(y-q), (7) 

and the goal of noise shaping amounts to arranging q in such a way that y — q is close to 
Ker(F). 

If the sequence (fj)™ of dual frame vectors were known to vary smoothly in j (including 
smooth termination into null vector), then EA quantization could be employed without 
much alteration, e.g., [5||2l]. However, this need not be the case for many examples of 
frames (together with their canonical duals) that are used in practice. For this reason, it 
has recently been proposed in to use special alternative dual frames, called Sobolev dual 
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frames, that are naturally adapted to £A quantization. It is shown in [4] (see also Section 
[2]) that for any frame E, if a standard rth order £A quantization algorithm with alphabet 
A = 5Z is used to compute q := g^A, then with an rth order Sobolev dual frame F := Fs b, r 
and xea := -Psobr^EA, the reconstruction error obeys the bound 



where D is the m x m difference matrix defined by 



1, if i = j, 

Dij := < -1, if i = j + l, (9) 
0, otherwise, 



and <7 m i n (.D r i?) stands for the smallest singular value of D r E. 
Contributions 

For the compressed sensing application that is the subject of this paper, E will simply be 
a sub-matrix of the measurement matrix <£>, hence it may have been found by sampling an 
i.i.d. random variable. Minimum singular values of random matrices with i.i.d. entries have 
been studied extensively in the mathematical literature. For anmxfc random matrix E with 
i.i.d. entries sampled from a sub-Gaussian distribution with zero mean and unit variance]^ 
one has 

a min (E) >^R-Vk (10) 
with high probability |26|. Note that in general D~ r E would not have i.i.d. entries. A naive 



lower bound for a min (D T E) would be a mm (D r )a min (E). However (see Proposition 3.1), 
Cmin(-D~ r ) satisfies 

(r miIl (D- r ) x r 1, (11) 

and therefore this naive product bound yields no improvement on the reconstruction error for 
SA-quantized measurements over the bound (HI) for PCM-quantized ones. In fact, the true 
behavior of cx mm (.D _r E) turns out to be drastically different and is described in Theorem |A| 



one of our main results (see also Theorem 3.7). 



For simplicity, we shall work with standard i.i.d. Gaussian variables for the entries of E. 
In analogy with our earlier notation, we define the "oversampling ratio" A of the frame E 
by 

A = p (12) 

Theorem A. Let E be an m x k random matrix whose entries are i.i.d. A/"(0, 1). For any 
a G (0, 1), if A > c{\ogm) l ^ l ~ a \ then with probability at least 1 — exp(— c'm\~ a ), 

a min (D- r E) > r A^)^, (13) 



As mentioned earlier, we do not normalize the measurement matrix $ in the quantization setting. 
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which yields the reconstruction error bound 

\\x-XxAh<r\- a{r -^S. (14) 

While the kind of decay in this error bound is familiar to £A modulation, the domain 
of applicability of this result is rather surprising. Previously, the only setting in which this 
type of approximation accuracy could be achieved (with or without Sobolev duals) was the 
case of highly structured frames (e.g. when the frame vectors are found by sampling along 
a piecewise smooth frame path). Theorem |A| shows that such an accuracy is obtained even 
when the analysis frame is a random Gaussian matrix, provided the reconstruction is done 
via Sobolev duals. 



In the compressed sensing setting, one needs (13) to be uniform for all the frames E that 
are found by selecting k columns of $ at a time. The proof of Theorem [X] extends in a 
straightforward manner using a standard "union bound" argument, provided A is known to 
be slightly larger. More precisely, if $ is an m x N matrix whose entries are i.i.d. according 



to W(0, 1), and if A := m/k > cQogiV) 1 ^ 1- "), then Q holds for all E = $ T with #T < k 
with the same type of probability bound (with new constants). This result can be utilized to 
improve the reconstruction accuracy of a sparse signal x from its EA-quantized compressed 
sensing measurements if the support T of x is known. This is because if T is known, $t is 
known, and its Sobolev dual can be found and used in the reconstruction. On the other hand, 
for most signals, recovering the exact or approximate support is already nearly guaranteed 
by the robust recovery result shown in ^ together with the stability of the associated SA 
quantizer. For example, a simple sufficient condition for full recovery of the support is that 
all the \xj\ for j G T be larger than C\\y — gsAlh for a suitable constant C. A precise version 
of this condition is stated in Theorem [Bl 

In light of all these results, we propose SA quantization as a more effective alternative 
of PCM (independent quantization) for compressed sensing. With high probability on the 
measurement matrix, a significant improvement of the reconstruction accuracy of sparse 
signals can be achieved through a two-stage recovery procedure: 

1. Coarse recovery: £i-minimization (or any other robust recovery procedure) applied 
to gsA yields an initial, "coarse" approximation of x, and in particular, the exact 
(or approximate) support T of x. 

2. Fine recovery: Sobolev dual of the frame $t applied to g^A yields a finer approxi- 
mation xsa of x. 



Combining all these, our second main theorem follows (also see Theorem 4.2): 

Theorem B. Let $ be an m x N matrix whose entries are i.i.d. according to jV(0, 1). 
Suppose a G (0,1) and A := m/k > c(log iV) 1 ^ 1-0 ^ where c = c(r,a). Then there are two 
constants d and C that depend only on r such that with probability at least 1 — exp(— c'm\~ a ) 
on the draw of $ 7 the following holds: For every x G such that minj esupp ^ \xj\ > C5, 
the reconstruction xea satisfies 

\\x-X S Ah<r \- a{r ~^5. (15) 
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To put this result in perspective, note that the approximation error given in ( 15 ) decays 
as the "redundancy" A = ? increases. In fact, by using an arbitrarily high order EA scheme, 
we can make this decay faster than any power law (albeit with higher constants). Note that 
such a decay is not observed in the reconstruction error bound for PCM given in (|4]). Of 
course, one could argue that these upper bounds may not reflect the actual behavior of 
the error. However, in the setting of frame quantization the performance of PCM is well 
investigated. In particular, let E be an m x k real matrix, and let K be a bounded set in 
For x G K, suppose we obtain <7pcm(#) by quantizing the entries of y — Ex using PCM 
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with alphabet A = b%. Let A opt be an optimal decoder. Then, Goyal et al. show in 17 
that 

[E||x-A opt (gp CM (x))||J 1/2 > \-H 

where A = m/k and the expectation is with respect a probability measure on x that is, 
for example, absolutely continuous. This lower bound limits the extent to which one can 
improve the reconstruction by means of alternative reconstruction algorithms from PCM- 
quantized compressed sensing measurements. On the other hand, setting, for example, 
a = 3/4 in Theorem [B] we observe that if we use a second-order EA scheme to quantize 
the measurements, and if we adopt the two-stage recovery procedure proposed above, the 
resulting approximation will be superior to that produced optimally from PCM-quantized 
measurements, provided m/k is sufficiently large. 

It is possible to imagine more sophisticated and more effective quantization and recov- 
ery algorithms for compressed sensing. However using EA quantization has a number of 
appealing features: 

• It produces more accurate approximations than any known quantization scheme in 
this setting (even when sophisticated recovery algorithms are employed). 

• It is modular in the sense that if the fine recovery stage is not available or practical 
to implement, then the standard (coarse) recovery procedure can still be applied as is. 

• It is progressive in the sense that if new measurements arrive (in any given order), 
noise shaping can be continued on these measurements as long as the state of the 
system (r real values for an rth order scheme) has been stored. 

• It is universal in the sense that it uses no information about the measurement matrix 
or the signal. 

The paper is organized as follows. We review the basics of EA quantization and Sobolev 
duals in frame theory in Section [2j followed by the reconstruction error bounds for random 
Gaussian frames in Section [3| We then present the specifics of our proposed quantization 
and recovery algorithm for compressed sensing in Section |4j We present our numerical 
experiments in Section [5] and conclude with extensions to more general settings in Section [6] 
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2 Background on HA quantization of frame expansions 



EA quantization 

The governing equation of a standard rth order £A quantization scheme with input y = (yj) 
and output q = (qj) is 

Vj q r j = 1,2,..., (16) 

where the qj G A are chosen according to some quantization rule given by 

Qj = Q( u j-i, • • • j u j-t, Vj, Vj-s)- (17) 

Not all EA quantization schemes are presented (or implemented) in this canonical form, 
but they all can be rewritten as such for an appropriate choice of r and u. We shall not 
be concerned with the specifics of the mapping Q, except that we need u to be bounded. 
The smaller the size of the alphabet A gets relative to r, the harder it is to guarantee this 
property. The extreme case is 1-bit quantization, i.e., = 2, which is typically the most 
challenging setting. We will not be working in this case. In fact, for our purposes, A will in 
general have to be sufficiently fine to allow for the recovery of the support of sparse signals. 
In order to avoid technical difficulties, we shall work with the infinite alphabet A = £Z, but 
also note that only a finite portion of this alphabet will be used for bounded signals. A 
standard quantization rule that has this "boundedness" property is given by the greedy rule 
which minimizes \uj\ given . . . ,Uj- r and yj, i.e., 

qj = argmm|^(-l) i ~ 1 (Mw i _ i + y i -a . (18) 

It is easy to check that with this rule, one has \uj\ < 2~ 1 5 and \yj — qj\ < 2 r ~ 1 5. In turn, 
if \\y\\oo < C, then one needs only L := 2\j~\ + 2 r + 1 levels. In this case, the associated 
quantizer is said to be log 2 L-bit, and we have 

\\u\\oo < 8 and \\y - qWoo < r 5. (19) 

With more stringent quantization rules, the first inequality would also have an r-dependent 
constant. In fact, it is known that for quantization rules with a 1-bit alphabet, this constant 



will be as large as 0(r r ), e.g., see 13,18 . In this paper, unless otherwise stated, we shall be 



working with the greedy quantization rule of ( 18 ) 



The initial condition of the recursion in ( 16 ) can be set arbitrarily, but it will be convenient 



for us to set them equal to zero for finite frames. With u_ r+ i = ■ ■ ■ = u = 0, and 



j = 1, . . . , m, the difference equation (16) can be rewritten as a matrix equation 



D r u = y — q, (20) 

where D is as in 

As before, we assume E is an m x k matrix whose rows form the analysis frame and 
F is a k x m left inverse of E whose columns form the dual (synthesis) frame. Given any 
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x G M. k , we set y = Ex, and define its rth order SA quantization g S A and its reconstruction 
£ea := Fqz&. Substituting (20) into ([7]), we obtain the error expression 

x-x = FD r u. (21) 

With this expression, \\x — x\\ can be bounded for any norm || • || simply as 

rn 

\\x-x\\ < H„£||(FI?0ill- (22) 
i=i 

Here (FD r )j is the jth column of FD r . This bound is also valid in infinite dimensions, 
and in fact has been used extensively in the mathematical treatment of oversampled A/D 
conversion of bandlimited functions. 

For r = 1, and the £2 norm, the sum term on the right hand side motivated the study of 
the so-called frame variation defined by 

m 

V(F):=J2Wh-fi+ih> (23) 
3=1 

where (fj) are the columns of F, and one defines f m +i — 0. Higher-order frame variations 
to be used with higher-order SA schemes are defined similarly, see |2}|3). Frames (analysis 
as well as synthesis) that are obtained via uniform sampling a smooth curve in M fc (so- 
called frame path) are typical in many settings. However, the "frame variation bound" is 
useful in finite dimensions when the frame path terminates smoothly. Otherwise, it does not 
provide higher-order reconstruction accuracy. Designing smoothly terminating frames can 
be technically challenging, e.g., |5). 

Sobolev duals 



Recently, a more straightforward approach was proposed in 21 for the design of (alternate) 
duals of finite frames for SA quantization. Here, one instead considers the operator norm 
of FD r on £2 and the corresponding bound 

\\x-x\\2<\\FD r \\ op \\u\\2. (24) 

Note that this bound is not available in the infinite dimensional setting of bandlimited 
functions due to the fact that u is typically not in £2. It is now natural to minimize ||FZ} r || op 
over all dual frames of a given analysis frame E. These frames, introduced in |3], have been 
called Sobolev duals, in analogy with £ 2 -type Sobolev (semi)norms. 

SA quantization algorithms are normally designed for analog circuit operation, so they 
control || ^ || 00, which would control ||u||2 only in a suboptimal way. However, it turns out 
that there are important advantages in working with the £2 norm in the analysis. The first 
advantage is that Sobolev duals are readily available by an explicit formula. The solution 
F so h, r of the optimization problem 

mm HFZrilop subject to FE = I (25) 
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is given by the matrix equation 



F soh ,.D r = {D- r E)\ (26) 

where ' stands for the Moore-Penrose inversion operator, which, in our case, is given by 
:= (E*E)~ 1 E* . Note that for r = (i.e., no noise-shaping, or PCM), one simply obtains 

F — E\ the canonical dual frame of E. 

The second advantage of this approach is that highly developed methods are present for 

spectral norms of matrices, especially in the random setting. Plugging (26) into (24), it 

immediately follows that 

\\x-x\\ 2 < ||(/r^)t|| op |H| 2 = J_ \\ u \\ 2 , (27) 

Vmm\L> E) 

where a m i n (D~ r E) stands for the smallest singular value of D~ r E. 



3 Reconstruction error bound for random frames 

In what follows, <Jj(A) will denote the jth largest singular value of the matrix A. Similarly, 
Xj{B) will denote the jth largest eigenvalue of the Hermitian matrix B. Hence, we have 
crj(A) = a/ Xj(A*A). We will also use the notation £(A) for the diagonal matrix of singular 
values of A, with the convention (E(A))jj = o~j(A). All matrices in our discussion will be 
real valued and the Hermitian conjugate reduces to the transpose. 

We have seen that the main object of interest for the reconstruction error bound is 
ffmm{D~ r E) for a random frame E. Let if be a square matrix. The first observation we make 
is that when E is i.i.d. Gaussian, the distribution of H(HE) is the same as the distribution 
of T,(Y,(H)E). To see this, let UH(H)V* be the singular value decomposition of H where U 
and V are unitary matrices. Then HE = UY,(H)V*E. Since the unitary transformation U 
does not alter singular values, we have H(HE) = T,(E(H)V*E), and because of the unitary 
invariance of the i.i.d. Gaussian measure, the matrix E := V*E has the same distribution as 
E, hence the claim. Therefore it suffices to study the singular values of H(H)E. In our case, 
H = D~ r and we first need information on the deterministic object S(Z)~ r ). The following 
result will be sufficient for our purposes: 

Proposition 3.1. Let r be any positive integer and D be as in @. There are positive 
numerical constants ci(r) and C2(r), independent of m, such that 

ci(r)(-J <o- 3 (D- r )<c 2 (r)(-) , j = l,...,m. (28) 

The proof of this result is rather standard in the study of Toeplitz matrices, and is given 
in Appendix [A} 
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3.1 Lower bound for <j m - m (D r E) 

In light of the above discussion, the distribution of a m ± n (D~ r E) is the same as that of 

inf \\£(D- r )Ex\\ 2 . (29) 

IM|2 = 1 

We replace H(D~ r ) with an arbitrary diagonal matrix S with Sjj =: Sj > 0. The first two 
results will concern upper bounds for the norm of independent but non-identically distributed 
Gaussian vectors. They are rather standard, but we include them for the definiteness of our 
discussion when they will be used later. 

Proposition 3.2. Let £ ~ A/"(0, ^I m ). For any 9 > 1, 

P (j>2 s)e 3 > ©INIL) < e m / 2 e-( @ - 1 ^ 2 . (30) 
Proof. Since Sj < ||s||oo for all j, we have 

(m \ / m \ 

E s % 2 > ii s i^J < p ^E^> e J- ( 31 ) 

This bound is the (standard) Gaussian measure of the complement of a sphere of radius 
\/mO and can be estimated very accurately. We use a simple approach via 

y>J>©) < min / e -(e-Er =1 ^ 2 )A/2TT e — ?/2 



o'=i / K j=i 

min e- A0/2 (l-A/m)- m / 2 

A>0 



= m/2 e- (e - 1)m/2 , (32) 
where in the last step we set A = m(l — -1 ). □ 

Lemma 3.3. Let E be an m x k random matrix whose entries are i.i.d. A/"(0, — ). For any 
> 1, consider the event 

S := |||S'-E'||^_ y gn < 2v / 0||s|| oo | . 

Then 

P(£ c ) < 5fce™/2 e -(e-i)m/ 2 _ 

Proof. We follow the same approach as in u\. The maximum number of p-distinguishable 
points on the unit sphere in M. k is at most (| + l) fc . (This follows by a volume argument^ 

c Balls with radii p/2 and centers at a p-distinguishable set of points on the unit sphere are mutually 
disjoint and are all contained in the ball of radius 1 + p/2 centered at the origin. Hence there can be at most 
(1 + p/2) k /(p/2) k of them. 
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as in e.g., 23, p. 487].) Fix a maximal set Q of —distinguishable points of the unit sphere 
in IR fc with j^Q < 5 k . Since Q is maximal, it is a |-net for the unit sphere. For each q G Q, 
consider = (Eq)j, j — 1, . . . , m. Then £ ~ A/"(0, — I m ). As before, we have 

m 

\\SEq\\l = J2^. 



3.2 



we have 



Let £(Q) be the event jdS-Eg^ < \/©||s||oo, Vg G Qj. Then, by Proposition 
the union bound 

P (£(QT) < 5 k Q m/2 e- {e ~ 1)m/2 . (33) 

Assume the event £(Q), and let M = \\SE\\ f k^ e m. For each ||x||2 = 1, there is q G Q with 
II? — x\\ 2 < 1/2, hence 

\\SEx\\ 2 < \\SEq\\ 2 + \\SE(x - q)\\ 2 < >/e||s|U + y ■ 

Taking the supremum over all x on the unit sphere, we obtain 

/— „ M 

M< ve||s||oo + y, 

i.e., HS'.E'll^^m < 2v^6||s|| 00 . Therefore C £, and the result follows. □ 



The following estimate concerns a lower bound for the Euclidean norm of (si£i, . . . , s m £ m ). 
It is not sharp when the Sj are identical, but it will be useful for our problem where Sj = 
(jj(D~ r ) obey a power law (see Corollary 3.5). 



Proposition 3.4. Let £ ~ A/"(0, ~Im)- For any 7 > 0, 

'e / ym\ L / 2 

3 3 I KL<m \ 

0=1 

Proof. For any t > and any integer L G {1, . . . , m}, we have 



' E«<^,s(t) ■.,)-. (34) 



I III \ n //t 



dx^ 



vi=i / i=i 

n 



a/27t Jm 



i 2 77 /" e -x2( m +t s 2) /2 dx 



< 



3=1 

e^n^^/m)- 1 / 2 
12 



< e t ^ 2 (m/t) L / 2 (s 1 s 2 ---s L )-\ (35) 

For any L, we can set t = L/7, which is the critical point of the function t t— > e ll t~ L . Since 

L is arbitrary, the result follows. □ 

Corollary 3.5. Let £ ~ A/"(0, — l m ), r be a positive integer, and c\ > be such that 



Sj > Ci 



m 
J . 



j = 1, . . . ,m. 



Then for any A > 1 and m > A, 

/ m 



(36) 



(37) 



VJ=1 



Proof. By rescaling s^-, we can assume C\ — \. For any L e {1, . . . , m}, we have 



(^•••^r 1 <^f<(8Lr/ 2/ Lr 



where we have used the coarse estimate L\ < e l/l2L (2nL) 1/2 (L/e) L < (8L) 1/2 (L/e) L . Setting 

r-1/2 



7 = A 2r 1 in Proposition 3.4 , we obtain 

/ m 

P S S ?^ <A2r_1 ) <( 8L ) r/2 



We set L = [_^J . Since 1 < A < m, it is guaranteed that 1 < L < m. Since AL < m, we get 

L 



AL\ 
em ) 



Plugging this in (38) and using 8e 2 < 60, we find 

p(X>|£-< A 21 - 1 ) < (eom/A)^-" 1 ^" 1 / 2 ^. 



(39) 



□ 

Theorem 3.6. Let E be an m ^ k random matrix whose entries are i.i.d. Af(Q, —), r be a 
positive integer, and assume that the entries Sj of the diagonal matrix S satisfy 

(40) 



ci y—j < Sj < c 2 m r , j = 1, . . . , m. 
Let A > 1 be any number and assume m > A. Consider the event 

T := ^\\SEx\\ 2 > ^ Cl A r - 1/2 ||x|| 2 , Vx e M fc j . 

Then 

P {JF C ) < 5 fc e" m / 2 + 8 r (17c 2 / Cl ) k A k ' 2 W r( * +1/2) e -m(r-l/2)A_ 
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Proof. Consider a p-net Q of the unit sphere of M. k with j^Q < (- + l) fe where the value 

of p < 1 will be chosen later. Let £{Q) be the event {||S£g|| 2 > ciA r " 1/2 , Vg G q}. By 
Corollary |3.5[ we know that 



P 



< 



9 



+ 1 



r/2 



-m(r-l/2)/A 



(41) 



3.3 



we 



Let £ be the event in Lemma 3.3 with = 4. Let E be any given matrix in th e ev ent 
£ n £(Q)- For each \\x\\2 = 1, there is g G Q with ||g — x\\% < p, hence by Lemma 
have 

|| ^(a; - g)|| 2 < 4HSIIOOH0; - g|| 2 < 4c 2 m r p. 



Choose 



Hence 



8c 2 m r 8c 2 VA^mJ 



c x /Ay 



1 



||^|| 2 > ||££g|| 2 - ||££(x - g)|| 2 > c^ 1 / 2 - 4c 2 m> = -CxA^. 

This shows that £ fl C J 7 . Clearly, p < 1/8 by our choice of parameters and hence 

- + 1 < Using the probability bounds of Lemma 3.3 and (41), we have 



P^) < 5 fc 4 m /V 3m/2 + ( 



'' / 60m^ ' 



A 



-m(r-l/2)/A 



< 5 V m / 2 + 8 r (17c 2 / Cl ) fc A fc / 2 (^) 



m x r(fc+l/2) 



-m(r-l/2)/A 



where we have used 2 < e and y60 < 8 for simplification. 

The following theorem is now a direct corollary of the above estimate. 



(42) 

□ 



Theorem 3.7. Let E be an m x k random matrix whose entries are i.i.d. A/"(0, — ) ; r be a 
positive integer, D be the difference matrix defined in and the constant c\ = Ci(r) be as 
in Proposition 3.1. Let < a < 1 be any number. Assume that 



ni 



A:=^>c 3 (logm) 1 /(i-«), 



where C3 = c%(r) is an appropriate constant. Then 

P (a min (/r r £) > ciA a(r '- 1/2) ) > 1 - 2e- Ciml ' ak ° 
for some constant C4 = c±(r) > 0. 



(43) 



(44) 
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Proof. Set A = X a in Lemma 3.6 We only need to show that 



max 



5 fc e- m/2 ,8 r (17c 2 / Cl ) fc A fc / 2 



m\ r (W/2) 



A 



-m(r-l/2)/A 



< e 



-C477X a fe° 



It suffices to show that 
and 



fc log 5 - m/2 < -c A m}- a k a 



1 1 L 777, 

rlog8 + fclog(17c 2 /ci) + -Hog A + r(fc + -) log(m/A) - ( r ~^)~^ - ~ c i m k a 



The first inequality is easily seen to hold if A > -j^sA. p or the second inequality, first notice 

2 C4 

that m/A = m l ~ a k a . Since fc + 1/2 x fc, and r — 1/2 x r, it is easily seen that we only need 
to check that 



k log m < c 5 



m 
A 



for a sufficiently small C5. This follows from our assumption on A by setting c$ = 1/c 



l-a 



□ 



Remark. By replacing i? in Theorem 3.7 with y/mE, we obtain Theorem [Aj 



3.2 Implication for compressed sensing matrices 



Theorem 3.8. Let r, D, cAr) be as in Theorem 3.7 and $ be an m x N random matrix 



whose entries are i.i.d. jV(0, -M. Lei < a < 1 be any number and assume that 

A:=^>c 6 (logA0 1/(1 - Q) , 



(45) 



where Cq = c§(r) is an appropriate constant. Then with probability at least 1 — 2e C7 ' mA " for 
some c-j = Cf{r) > 0, every m x k submatrix E of satisfies 



a miQ (D- r E) > aX^-W. 



(46) 



Proof. We will choose C7 = C4/2, where C4 is as in Theorem 3.7 The proof will follow 
immediately by a union bound once we show that 



N 



Since ( fc ) < N k , it suffices to show that 



klogN < -m l - a k a . 
& - 2 



Both this condition and the hypothesis of Theorem 3/7 will be satisfied if we choose 

c 6 = max( C3 ,(2/c 4 ) 1/{1 - Q) ). 



□ 
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Remark. If $ is a Gaussian matrix with entries i.i.d. W(0, 1) rather than jV(0, — ), Theo- 
applied to ^=< 



rem 



3.8 



1 =$ implies that every m x k submatrix E of $ satisfies 



a min (D- r E) > Cl A^ 



-1/2) 



m. 



(47) 



4 EA quantization of compressed sensing measurements 

In this section we will assume that the conditions of Theorem 13.81 are satisfied for some 



< a < 1 and r, and the measurement matrix $ that is drawn from Af(0, 1) yields (47). 
For definiteness, we also assume that $ admits the robust recovery constant C\ = 10, i.e., 
the solution x # of the program ([I]) satisfies 



\y — y\\2<e =>• \\x - x*\\ 2 < 10— !=e. 



Note again that our choice of normalization for the measurement matrix $ is different 
from the compressed sensing convention. As mentioned in the Introduction, it is more 
appropriate to work with a measurement matrix $ ~ A/"(0, 1) in order to be able to use a 
quantizer alphabet that does not depend on m. For this reason, in the remainder of the 
paper, $ shall denote an m x N matrix whose entries are i.i.d. from A/"(0, 1). 

Let q := g^A be output of the standard greedy rth order EA quantizer with the alphabet 
A = blj and input y. As stated in Section [2j we know that \\y — qW^ < 2 r ~ 1 5 and therefore 
\\y — <i\\2 < 2 r ~ l 5\fm. 



Coarse recovery and recovery of support 

Our first goal is to recover the support T of x. For this purpose we shall use a coarse 
approximation of x. Let 

x := argmin \\z\\\ subject to ||<&2 — q\\ 2 < e := 2 r ^ 1 Sy/m. (48) 

By the robust recovery result (for our choice of normalization for $), we know that 

\\x — x'\\2 < J) := 5 • 2 r S. 

The simplest attempt to recover T from x' is to pick the positions of its k largest entries. 
This attempt can fail if some entry of Xj on T is smaller than rj for then it is possible that 
x'j = and therefore j is not picked. On the other hand, it is easy to see that if the smallest 
nonzero entry of x is strictly bigger than 2rj in magnitude, then this method always succeeds. 
(Since \\x — x'W^ < rj, the entries of x' are bigger than rj on T and less than rj on T c .) The 
constant 2 can be replaced with \/2 by a more careful analysis, and can be pushed arbitrarily 
close to 1 by picking more than k positions. The proposition below gives a precise condition 
on how well this can be done. We also provide a bound on how much of x can potentially 
be missed if no lower bound on \xj\ is available for j G T. 
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Proposition 4.1. Let \\x — x'\\ £ n < r\, T — supp x and k = \T\. For any k' G {k, . . . ,N- 

let T' be the support of (any of) the k' largest entries of x' . 

(i) \\x T \ T ,\\ 2 < /3rj where /3 < (l + ^) 1/2 . 

(ii) If \xj\ > 777 for all j G T, where 7 := (l + fc ,_ 1 fc+1 ) 1 ^ 2 , then T' D T. 
Proof, (i) We have 

El / 1 2 1 V"^ I / |2 11 / 1 1 2 2 

\ X j ~ X j\ + y, \ X j\ = \\ x — x \\2 — V ■ 

In particular, this implies 



j£T\T' jST'\T 



2 ^ 2 



Suppose T \ T" 7^ 0. Then X" \ T is also nonempty. In fact, we have 

|T"\T| = \T\T'\ + k' - k. 

Now, observe that 

1 \ ^ I / |2 ^ I / |2 • I /■ 1 2 1 \ ^ I / i2 

r— r-— -; > \xA < max a;,- < mm \xA < lrTn . — -, > , 
\T\T'\ ^ 1 Jl ~ jer\T> 1 Jl _ ier'\T 1 jl ~ \T'\T\ 1 jl 



which, together with ( 50 ) implies 



||z r \T'||2 < lkr\T'H 2+ W( X ~ X ')T\T'h < \\ x T\T'h + \/v 2 ~ |y ^| lkr\ T /IH- 



It is easy to check that for any A > 0, and any < t < rj/y/ A, 

1/2 



t+y/ V *- At 2 < 



The result follows by setting A = |T' \ T|/|T \ T'\ and noticing that A > k'/k. 

(ii) Let zi > ■ ■ ■ > zjq be the decreasing rearrangement of |a^|, . . . , |a^y|. We have 



jeT i=i 



z. 2 



so 

AT fe'+l 



EKi 2 ^ E ^> E ^>(k'-k+i)4 +1 . 

j£T c i=k+l i=k+l 
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Hence by ( 49 ) we have 



max \xj — x'j\ 2 + (k' — k + l)z k , +1 < r/ z 
Since \x'j\ > \ Xj I I X j X a I, the above inequality now implies 



min \x'A > min \xA — max \x~ — x'A > min \xA — \/n 2 — (k' — k + \)zh<-\- 
jeT 1 3 1 ~~ jeT 1 3 1 ier 1 J 3 1 ~~ jeT 171 V v ' k +i 



Now, another application of (51 ) with A = k! — k + 1 yields 



and therefore 

min I a;'. I > min |se,-| + Zy+i — jrt > zy+\ = max \xA 
jeT 3 jer J jeT'c 3 



It is then clear that TcT' because if T' c fl T ^ 0, the inequality 

max I x'- 1 > max Ix'-I >min|x'-| 
jeT' c 3 jeT' c nT 7 jeT J 

would give us a contradiction. □ 

Note that if the k' largest entries of x' are picked with k' > k, then one would need to 
work with X" for the fine recovery stage, and therefore the starting assumptions on $ have 
to be modified for k' . For simplicity we shall stick to k' — k and consequently 7 = \/2. 



Fine recovery 

Once T is found, the rth order Sobolev dual frame F := -Fsob,r of E = $^ is computed and 
we set xsa = Fq. We now restate and prove Theorem [Bj 

Theorem 4.2. Let $ be an m x N matrix whose entries are i.i.d. according to jV(0, 1). 
Suppose a £ (0,1) and A := m/k > c(log jV) 1 ^ 1-0 ^ where c = c(r, a). Then there are two 
constants d and C that depend only on r such that with probability at least 1 — exp(— c'm\~ a ) 
on the draw of $ ; the following holds: For every x G such that mhijgsupp^) \xj\ > C8, 
the reconstruction xea satisfies 

\\x-x^\\ 2 < r \- a ^5. (52) 



Proof. Suppose that A > c(log N) 1 ^ 1 °^ with c = cq as in the proof of Theorem 3.8 Let gsA 



be obtained by quantizing y := via an rth order SA scheme with alphabet A = 5Z, and 



with the quantization rule as in (18), and let u be the associated state sequence as in (16) 



Define as the solution of the program 

min ||z||i subject to \\<&z — gsA||2 < £■ 
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Suppose that $ admits the robust recovery constant Ci, i.e., the solution of the program 
(|3j) satisfies \\x — x # ||2 < Cie/ -^/m for every x in provided that \\y — <7ea|| < £• Note that 
Ci, as given for example in [9], only depends on the RIP constants of $ and is well-behaved 
if m and A" satisfy the hypothesis of the theorem. As discussed in Section [2j in this case we 
have II y — gsAlh < 2 r ~ 1 5 A /m which implies 



||a;-a; # || 2 < d^S. 

Assume that 



mm \ x 

ier 



j\ >C l -2 r - 1 / 2 5=:C5. (53) 



Then, Proposition 4.1 (with 7 = \/2 and 77 = C\l r 1 ) shows that T', the support of the k 



largest entries of x # , is identical to the support T of x. Finally, set 

^EA = -Psob.r^SA 

where -F sobr is the rth order Sobolev dual of $t- Using the fact that ||w|| 2 < 2^ l 5y/m 
(see Section [2]) together with the conclusion of Theorem 3.8 and the error bound (27), we 
conclude that 

Hulls A -a(r ~ 1/2) 
y/ma min {D r E) 2ci 

Note that the RIP and therefore the robust recovery will hold with probability 1 — exp(c"m), 
and our Sobolev dual reconstruction error bound will hold with probability 1— exp(— crm\^ a ). 



Here c\ and Cj are as in the proof of Theorem 3.8 



□ 

Remark. To interpret the size condition in a concrete case, assume that $ admits the robust 
recovery constant C\ = 10, and that we have 

min \xA > V2n = 5 ■ 2 r+1/2 5. (55) 

If PCM is used as the quantization method, then the best error guarantee we have that holds 
uniformly on T would be 

\\x -^pcmII 00 — 11"^ ^pcmII 2 — 55. 

It can be argued that the approximately recovered entries of Xp CM are meaningful only when 
the minimum nonzero entry of x is at least as large as the maximum uncertainty in Xp CM , 



which is only known to be bounded by 55. Hence, in some sense the size condition (55) is 
natural (modulo the factor 2 r+1//2 ). 
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Quantizer choice and rate-distortion issues 

So far we have not made any assumptions on the step size 5 of the uniform infinite quantizer 
A = 5Z. An important question concerns how large S should be for the most effective use of 
resources. This question is motivated by the fact that infinite quantizers are not practical 
and have to be replaced by finite ones. In the same vein, an alternative question is to 
determine the minimum number of bits that the quantizer needs to incorporate as well as 
the resulting approximation error. First, let us assume that 

A < \xj\ < p := 2 b A for all j e T. (56) 

For usefulness of our results, one would be interested in the regime A <C /). Thus, we 
introduce 2 b = p/A to represent the number of dyadic scales over which the input is allowed 
to range. Clearly, 8 r , the quantization step size used by an rth order SA scheme for our 
support recovery results to hold must satisfy 5 r < 2 f+i /2 (as before, we assume C\ = 10). 
Let us for the moment use the largest allowable step-size, i.e., set 

' - A " (57) 



2V+1/2- 



Next, let us assume that a B r -bit uniform quantizer of step size 5 r is to replace A = <5Z. 
We know that ||g||oo could be as large as 2 r ~ 1 5 r + ||y||oo, therefore we need to bound \\y\\oo 
efficiently If we use the RIP, then $ does not expand the ^ 2 -norm of /c-sparse vectors by 
more than a factor of 2y / m (note our choice of normalization for $), and therefore it follows 
that 

||y||oo < IMI2 < 2v / m||a;||2 < 2pVmk, 
which is a restatement of the inequality 



\E\U^ < y/k\\E\ 



2 



that holds for any m x k matrix E. However, it can be argued that the (00, oo)-norm of 
a random matrix should typically be smaller. In fact, if E were drawn from the Bernoulli 
model, i.e., Eij ~ ±1, then we would have 

\\E\\fk .cm = k = A~ 1//2 v / mA;, 
as can easily be seen from the general formula 

k 

KKra ' ' 
3=1 

Using simple concentration inequalities for Gaussian random variables, it turns out that 
for the range of aspect ratio A = m/k and probability of encountering a matrix $ that we 
are interested in, we have \\E\\ < \- a / 2 y/mk for every m x k submatrix E of $. We 

start with the following estimate: 
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Proposition 4.3. Let £i> ■ ■ ■ >£fc standard Gaussian variables. Then, for any > 1, 



pj^y >e) < 2V 



(59) 



Proof. 



p(E^i >e ) ^ 



mm 

t>o 



k 

f e -(e-E* =1 l^l)tTT e - 



2/2 dXj 



2tt 



mm e 

t>o 



dx 
s/2tt 



mm 

t>o 



in e- e M2e* / 2 / e"* 



~ '> l°° -h(x-t) 2 dx 



>2tx 



< 2 fc min e - @t+kt2 ' 2 
t>o 

= 2 fc e _6)2/{2fc) 



where in the last step we set t = Q/k. 



(60) 

□ 



Proposition 4.4. Let $ be an m x N random matrix whose entries are i.i.d. jV(0,1). Let 
< a < 1 fre any number and assume that 



m 



\:=j>*{logN) 1 '<- 1 - a \ 



(61) 



where C\ is an appropriate constant. Then with probability at least 1 — e C2ml afeQ for some 
c 2 > 0, every m x k submatrix E of $ satisfies 



\E\ 



Proof. Proposition 4.3 straightforwardly implies that 



P ({3T such that \T\ = k and \\^ T \\^^ > ©}) < KJ m2 k e'® 2 ' {2k) . 



(62) 



(63) 



Let = X a / 2 \/ mk. It remains to show that 



k log N + k log 2 + log m + c 2 m L ~ a k a < 



2Jfe" 



If Ci in (61) is sufficiently large and C2 is sufficiently small, then the expression on the left 
hand side is bounded by k\ 1 ~ a /2 = Q 2 /(2k). □ 

Without loss of generality, we may now assume that $ also satisfies the conclusion of 



Proposition |4.4| Hence we have an improved bound on the range of y given by 

IMIoc < p\~ a/2 V^n~k = pA {1 - Q)/2 fc. (64) 
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We assume B r is chosen to satisfy 



)B r -l 



6 r = 2 r ~ 1 5 r + pA (1 " Q)/2 A;, 



(65) 



so that the quantizer is not overloaded. Since p/S r ~ 2 r+l l 2+b by (56) and (57), we see that 
the second term on the right hand side of (65) is significantly larger than the first, which 
implies 

2 B '- 1 5 r « 2 b AX {1 ~ a)/2 k. (66) 



Hence, using (57) again, B r must satisfy 

2-Br-l pj, 5 2 6+r+l/2^(l-a)/2^, 



(67) 



Based on Theorem 4.2, the approximation error (the distortion) T)-^a incurred after the fine 
recovery stage via Sobolev duals satisfies the bound 



£SA <r \- a(r - 1/2) 5 r 



2 r + 1 /2 ' 



(68) 



A similar calculation for the PCM encoder with the same step size S r and the standard 
i\ decoder results in the necessity for roughly the same number of bits B r as the E A encoder 
(because of the approximation (66)), but provides only the distortion bound 

.4 



l>CM < S r ~ 



2 r+l/2 ' 



(69) 



Note that the analysis above requires that both PCM and £A encoders utilize high- 
resolution quantizers, however the benefit of using EA encoders is obvious upon comparing 
([68l and (ph. 



5 Numerical experiments 

In order to test the accuracy of Theorem |3.7[ our first numerical experiment concerns the 
minimum singular value of D~ r E as a function of A = m/k. In Figure [TJ we plot the worst 
case (the largest) value, among 1000 realizations, of l/a m i n (D~ r E) for the range 1 < A < 25, 
where we have kept k = 50. As predicted by this theorem, we find that the negative slope 
in the log- log scale is roughly equal to r — 1/2, albeit slightly less, which seems in agreement 
with the presence of our control parameter a. As for the size of the r-dependent constants, 
the function 5 r A~ r+1 / 2 seems to be a reasonably close numerical fit, which also explains why 
we observe the separation of the individual curves after A > 5. 

Our next experiment involves the full quantization algorithm for compressed sensing in- 
cluding the "recovery of support" and "fine recovery" stages. To that end, we first generate 
a 1000 x 2000 matrix 4>, where the entries of $ are drawn i.i.d. according to A/"(0, 1). To 
examine the performance of the proposed scheme as the redundancy A increases in compar- 
ison to the performance of the standard PCM quantization, we run a set of experiments: In 
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each experiment we fix the sparsity k G {5, 10,20,40}, and we generate fc-sparse signals x 
with the non-zero entries of each signal supported on a random set T, but with magnitude 
1/y/k. This ensures that ||x|| 2 = 1. Next, for m G {100, 200, 1000} we generate the mea- 
surements y = Q^x, where is comprised of the first m rows of $. We then quantize y 



using PCM, as well as the 1st and 2nd order SA quantizers, defined via (16) and (18) (in all 
cases the quantizer step size is 5 = 10 -2 ). For each of these quantized measurements q, we 
perform the coarse recovery stage, i.e., we solve the associated l\ minimization problem to 
recover a coarse estimate of x as well as an estimate T of the support T. The approximation 
error obtained using the coarse estimate (with PCM quantization) is displayed in Figures [2] 
and [3] (see the dotted curve). Next, we implement the fine recovery stage of our algorithm. 
In particular, we use the estimated support set T and generate the associated dual F sob r . 
Defining -F so h>n : = 

($^ m) )t, in each case, our final estimate of the signal is obtained via the 
fine recovery stage as Xf = F so ^ r q, Xf c = 0. Note that this way, we obtain an alternative 
reconstruction also in the case of PCM. We repeat this experiment 100 times for each (k, m) 
pair and plot the average of the resulting errors \\x — a: || 2 as a function of A in Figure [2] as well 
as the maximum of ||x — ^ II 2 m Figure [3] For our final experiment, we choose the entries of 
xt i.i.d. from A/"(0, 1), and use a quantizer step size 5 = 10~ 4 . Otherwise, the experimental 
setup is identical to the previous one. The average of the resulting errors \\x — x\\2 as a 
function of A is reported in Figure [4] and the maximum of \\x — x\\2 in Figure [5] 
The main observations that we obtain from these experiments are as follows: 

• SA schemes outperform the coarse reconstruction obtained from PCM quantized mea- 
surements significantly even when r = 1 and even for small values of A. 

• For the SA reconstruction error, the negative slope in the log-log scale is roughly equal 
to r. This outperforms the (best case) predictions of Theorem [B] which are obtained 
through the operator norm bound and suggests the presence of further cancellation 
due to the statistical nature of the SA state variable u, similar to the white noise 
hypothesis. 

• When a fine recovery stage is employed in the case of PCM (using the Moore-Penrose 
pseudoinverse of the submatrix of $ that corresponds to the estimated support of x), 
the approximation is consistently improved (when compared to the coarse recovery). 
Moreover, the associated approximation error is observed to be of order 0(A -1 / 2 ), 
in contrast with the error corresponding to the coarse recovery from PCM quantized 
measurements (with the l\ decoder only) where the approximation error does not 
seem to depend on A. A rigorous analysis of this behaviour will be given in a separate 
manuscript. 
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6 Remarks on extensions 



6.1 Other noise shaping matrices 

In the above approach, the particular quantization scheme that we use can be identified 
with its "noise-shaping matrix" , which is D r in the case of an rth order E A scheme and the 
identity matrix in the case of PCM. 

The results we obtained above are valid for the aforementioned noise-shaping matrices. 
However, our techniques are fairly general and our estimates can be modified to investigate 
the accuracy obtained using an arbitrary quantization scheme with the associated invertible 
noise-shaping matrix H. In particular, the estimates depend solely on the distribution of the 
singular values of H. Of course, in this case, we also need change our "fine recovery" stage 
and use the "if-dual" of the corresponding frame E, which we define via 



As an example, consider an rth order high-pass EA scheme whose noise shaping matrix 
is H r where H is defined via 



It is easy to check that the singular values of H are identical to those of D. It follows that 
all the results presented in this paper are valid also if the compressed measurements are 
quantized via an an rth order high-pass EA scheme, provided the reconstruction is done 
using the if r -duals instead of the rth order Sobolev duals. Note that such a result for 
high-pass EA schemes is not known to hold in the case of structured frames. 

6.2 Measurement noise and compressible signals 

One of the natural questions is whether the quantization methods developed in this paper 
are effective in the presence of measurement noise in addition to the error introduced during 
the quantization process. Another natural question is how to extend this theory to include 
the case when the underlying signals are not necessarily strictly sparse, but nevertheless still 
"compressible" . 

Suppose x G WL N is not sparse, but compressible in the usual sense (e.g. as in |9j), and let 
y = $x+e, where e stands for additive measurement noise. The coarse recovery stage inherits 
the stability and robustness properties of £% decoding for compressed sensing, therefore the 
accuracy of this first reconstruction depends on the best fc-term approximation error for x, 
and the deviation of $x from the quantized signal q (which comprises of the measurement 
noise e and the quantization error y — q). Up to constant factors, the quantization error for 
any (stable) EA quantizer is comparable to that of PCM, hence the reconstruction error at 
the coarse recovery stage would also be comparable. In the fine recovery stage, however, the 
difference between c m£lx (F H H) and cr max (F H ) plays a critical role. In the particular case of 



F H H = (HE)*. 



(70) 




1, if i = j or if i = j + 1, 
0, otherwise. 



(71) 
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H = D r and F H = F sob r , the Sobolev duals we use in the reconstruction are tailored to reduce 
the effect of the quantization error introduced by an rth order EA quantizer. This is reflected 
in the fact that as A increases, the kernel of the reconstruction operator -F so b,r contains a 
larger portion of high-pass sequences (like the quantization error of SA modulation), and is 
quantified by the bound a max (F soh , r D r ) < X-^l^mT 1 ^ (see Theorem |X| @ and (|27}). 
Consequently, obtaining more measurements increases A, and even though | y — q\\ 2 increases 
as well, the reconstruction error due to quantization decreases. At the same time, obtaining 
more measurements would also increase the size of the external noise e, as well as the "aliasing 
error" that is the result of the "off-support" entries of x. However, this noise+error term 
is not counteracted by the action of F so ^ r . In fact, for any dual F, the relation FE = I 
implies cr max (F) > l/cr max (E) > vrC 1 ! 2 already and in the case of measurement noise, it is 
not possible to do better than the canonical dual E^ on average. In this case, depending on 
the size of the noise term, the fine recovery stage may not improve the total reconstruction 
error even though the "quantizer error" is still reduced. 

One possible remedy for this problem is to construct alternative quantization schemes 
with associated noise-shaping matrices that balance the above discussed trade-off between 
the quantization error and the error that is introduced by other factors. This is a delicate 
procedure, and it will be investigated thoroughly in future work. However, a first such 
construction can be made by using "leaky" SA schemes with H given by 

{1, if * = j, 
-fi i£i = j + l, (72) 
0, otherwise, 

where fj, E (0, 1). Our preliminary numerical experiments (see Figures [6] and [7j suggest that 
this approach can be used to improve the accuracy of the approximation further in the fine 
recovery stage in this more general setting. We note that the parameter \i above can be 
adjusted based on how compressible the signals of interest are and what the expected noise 
level is. 
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A Singular values of D r 

It will be more convenient to work with the singular values of D T '. Note that because of our 
convention of descending ordering of singular values, we have 

^{D- r ) = 1 -jj^- v J = l,...,m. (73) 

For r = 1, an explicit formula is available [27l|28]. Indeed, we have 



which implies 



0j(D- 1 ) = 7 !■ J = l,-..,m. (75) 



2 S in ( 7r ^~ 1 / 2 ) ^ 

^ b111 \2(m+l/2) J 



The first observation is that o~j(D r ) and (<jj(D)) r are different, because .D and D* do 
not commute. However, this becomes insignificant as m — > oo. In fact, the asymptotic 
distribution of (<yj(D r ))™ =1 as m — » oo is rather easy to find using standard results in the 
theory of Toeplitz matrices: D is a banded Toeplitz matrix whose symbol is f(9) = 1 — e* 6 *, 
hence the symbol of D r is (1 — e ie ) r . It then follows by Parter's extension of Szego's theorem 



25 that for any continuous function ip, we have 

hm -f>fo(£r)) = ^ r nf{9)\ r )de. (76) 

3=1 7_7r 

We have |/(6 I )| = 2sin|^|/2 for \9\ < it, hence the distribution of (<7 7 (-D r ))™ =1 is asymp- 
totically the same as that of 2 r sin r (7rj/2m), and consequently, we can think of aj(D~ r ) 
roughly as (2 r sin r (7T?'/2m)) . Moreover, we know that ||D r || op < ||-D||o P < 2 r , hence 
cr miQ (D- r )>2- r . 

When combined with known results on the rate of convergence to the limiting distribution 
in Szego's theorem, the above asymptotics could be turned into an estimate of the kind given 
in Proposition 3.1 perhaps with some loss of precision. Here we shall provide a more direct 
approach which is not asymptotic, and works for all m > 4r. The underlying observation is 
that D and D* almost commute: D*D — DD* has only two nonzero entries, at (1, 1) and 
(m,m). Based on this observation, we show below that D* r D r is then a perturbation of 
(D*D) T of rank at most 2r. 

Proposition A.l. Let = D* r D r — (D*D) r where we assume m > 2r. Define 
I r := {l,...,r} x {l,...,r}U {m — r + 1, . . . , m} x {m — r + 1, . . . , m}. 
Then CffJ = for all G 1%. Therefore, rank(C^) < 2r. 
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Proof. Define the set G r of all "r-cornered" matrices as 

e r = {M:M M = 0if {i,j)el?}, 
and the set 23 r of all "r-banded" matrices as 

■B r = {M : M id = if \i - j\ > r}. 

Both sets are closed under matrix addition. It is also easy to check the following facts (for 
the admissible range of values for r and s): 

(i) If B E 23 r and C E G s , then BC E G r+S and CB E G r+S . 

(ii) If B E !B r and B E S s , then BB E T> r+S . 

(iii) If C E G r and C E G s , then CC E e max ( rjS ). 

(iv) If C E C r , then D*CD E G r+1 . 

Note that DD*,D*D E 'B 1 and the commutator [D*,D] =: T 1 E Ci. Define 

T r := (D*D) r - (DD*) r = (DD* + - (DD*) r . 

We expand out the first term (noting the non-commutativity), cancel (DD*) r and see that 
every term that remains is a product of r terms (counting each DD* as one term) each of 
which is either in 23 1 or in Ci. Repeated applications of (i), (ii), and (iii) yield T r E G r . 

We will now show by induction on r that C^ r > E G r for all r such that 2r < m. The cases 
r = and r = 1 hold trivially. Assume the statement holds for a given value of r. Since 

C (r+i) = D*(C {r) + T r )D 
and T r E G r , property (iv) above now shows that C^ r+ ^ E G r+1 . □ 



The next result, originally due to Weyl (see, e.g., 19, Thm 4.3.6]), will now allow us to 
estimate the eigenvalues of D* r D r using the eigenvalues of (D*D) r : 

Theorem A. 2 (Weyl). Let B and C be m x m Hermitian matrices where C has rank at 
most p. Then 

X j+P (B) < \{B + C)< \j- p (B), j=p + l,...,m-p, (77) 
where we assume eigenvalues are in descending order. 



We are now fully equipped to prove Proposition 3.1 
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Proof of Proposition^^ We set p = 2r, B = (D*D) r , and C = C (r) = D* r D r - (D*D) r in 
Weyl's theorem. By Proposition |A.l , C has rank at most 2r. Hence, we have the relation 

X j+2r ((D*D) r ) < Xj(D* r D r ) < A i _ 2r ((£>*£>) r ), j = 2r + 1, . . . , m - 2r. (78) 

Since Xj((D*D) r ) = Xj(D*D) r , this corresponds to 

aj + 2r(D) r < o-j(D r ) < o-j_2r(-D) r , j = 2r + 1, . . . , m - 2r. (79) 

For the remaining values of j, we will simply use the largest and smallest singular values of 
D r as upper and lower bounds. However, note that 



and similarly 



\D r \\ < \\D\\ r 

\u ||op _i ll-^llop 



lollop > \\D 



-1 II -r 

Hop 



Hence ( 79 ) can be rewritten as 

&min(j+2r,m)(D) r < 0~j(D T ) < CT m ax(j-2r,l) 



°i{D)) r 
(D) r } j = 1, 



, m. 



Inverting these relations via (|73j), we obtain 

) (D- 1 ) r < a j (D- r ) < a 



a 



mm(j+2r,m) 



max(j— 2r,l) 



-l\r 



1, m. 



(80) 



Finally, to demonstrate the desired bounds of Proposition 3.1, we rewrite (75) via the in- 
equality 2x/n < sinx < x for < x < tt/2 as 



m + 1 / 2 ^/n-iw ^ + 1/2 



(82) 



tt(; -1/2) " JV / -2(j-l/2)' 

and observe that min(j + 2r, m) x r j and max(j — 2r, 1) x r j for j = 1, . . . , m. □ 

Remark. The constants Ci(r) and Oz{r) that one obtains from the above argument would 
be significantly exaggerated. This is primarily due to the fact that Proposition 3.1 is not 
stated in the tightest possible form. The advantage of this form is the simplicity of the 
subsequent analysis in Section 3.1 Our estimates of a n { m {D~ r E) would become significantly 
more accurate if the asymptotic distribution of aj(D~ r ) is incorporated into our proofs in 
Section 3.1 However, the main disadvantage would be that the estimates would then hold 
only for all sufficiently large m. 
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X = m/k 



Figure 1: Numerical behavior (in log-log scale) of l/a m i n (D~ r E) as a function of A = m/k, 
for r = 0, 1, 2, 3, 4. In this figure, k = 50 and 1 < A < 25. For each problem size, the largest 
value of 1 / a m i n (D~ r E) among 1000 realizations of a random m x k matrix E sampled from 
the Gaussian ensemble 7V(0, ^.I m ) was recorded. 
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performance of various quantization/decoding schemes, k = 5 



performance of various quanf izalion/decoding schemes, k = 1 
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Figure 2: The average performance of the proposed EA quantization and reconstruction 
schemes for various values of k. For this experiment the non-zero entries of x are constant 
and 5 = 0.01. 
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performance of various quantization/decoding schemes, k = 5 



performance of various quanf izalion/decoding schemes, k = 1 
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Figure 3: The worst case performance of the proposed SA quantization and reconstruction 
schemes for various values of k. For this experiment the non-zero entries of x are constant 
and 5 = 0.01. 
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performance of various quantization/decoding schemes, k = 5 



performance of various quantization/decoding schemes, k = 10 
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Figure 4: The average performance of the proposed SA quantization and reconstruction 
schemes for various values of k. For this experiment the non-zero entries of x are i.i.d. 
JV(0, 1) and 5 = 10" 4 . 
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Figure 5: The worst case performance of the proposed SA quantization and reconstruction 
schemes for various values of k. For this experiment the non-zero entries of x are i.i.d. 
JV(0, 1) and 5 = 10" 4 . 
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Figure 6: The average case performance of the proposed SA quantization and reconstruction 
schemes (with general duals) for various values of k. For this experiment the non-zero entries 
of x are constant and 5 = 0.01. 
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performance of various quantization/decoding schemes, k = 5 



performance of various quaniization/decoding schemes, k = 1 
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Figure 7: The worst case performance of the proposed SA quantization and reconstruction 
schemes (with general duals) for various values of k. For this experiment the non-zero entries 
of x are constant and 5 = 0.01. 
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